% inputs:
% gm: R*p matrix of weights, each row representing different values of gamma
% u : n*K matrix of generalized residuals
% w : n*p matrix of instruments
% lambda : scalar penalty level

% outputs:
% ftn    : the value of the test statistic

function [ftn, output] = hdm(gm, u, w, lambda, varargin)
ip = inputParser;
addParameter(ip, 'dm', 0);
parse(ip, varargin{:});
dm = ip.Results.dm;  % indicator for exponential weight demeaning.

p = size(gm, 2);
n = size(u, 1);
K = size(u, 2);

if isequal(dm, 1)
    mw = mean(exp(w * gm'));
    W = (exp(w * gm') - mw);
    num = u' * W;    % K by R matrix
                     % Difference from "hdm".
    den = (u.^2)' * (W.^2);
else
    num = (u)' * (exp(w * gm'));
    den = (u.^2)' * (exp(w * gm')).^2;
end

output.Q = - abs(num) ./ sqrt(den);
output.l1 = sum(abs(gm'), 1);
ftn = output.Q + lambda * repmat(output.l1, K, 1);
end